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^C) Abstract 

^H 

Boltzmann samplers introduced by Duchon et al. in 2001 is a method by which to automatically 

£N) build random samplers from a combinatorial class' symbolic specification. But after over a decade, no 

complete standalone implementation of these samplers has been achieved. This paper describes how 

Qh to address the primary roadbloack: that the concept at the heart of Boltzmann sampling, evaluating 

generating functions which are sometimes only known implicitly, is difficult to implement automatically 

v<-s and efficiently enough. 

By adapting the rejection method, a classical tool from the random sampling literature, we show how 

to bypass this issue almost entirely. We then redefine the Boltzmann sampler framework and show how to 

build concretely useful and efficient algorithms. Our work radically simplifies much of the computational 

burden that was previously associated with Boltzmann samplers, and hints at future powerful extensions. 

CO 

i — i 1 Introduction 



Random generation of large combinatorial objects is a fundamental problem with countless applications 
in scientific modeling. This paper, for the first time, outlines the basis for the development of a complete 
system that can automatically solve this problem for objects comprised of millions or billions of atoms. 

The topic of this paper is the efficient and automatic construction of random samplers which uniformly 
draw large combinatorial objects from among objects of the same size. 

Striving for uniformity is unintuitive and difficult. The only combinatorial object that is naturally easy 
to generate uniformly are permutations. 
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Automatic combinatorial samplers. Nijenhuis and Wilf introduced the recursive method [12] in the 
late 70s (later extended by [9]), the first automatic random generation method; it is termed automatic be- 
cause it can directly derive random samplers from any combinatorial description — no bijection, no clever 
algorithm, no complicated equations are needed. The drawback is that this method is costly: to generate an 
object of size n, it requires knowing the complete enumeration of the combinatorial class up to size n; and 
predictably when n is large, this enumeration is significant both to calculate and to store. 

Enter Boltzmann sampling, introduced by Duchon et al. in 2002 O [5), of which the key insight was 
that the coefficients do not need to be extracted: instead, correct probabilities can be obtained by proxy, by 
evaluating the counting generating functions, that is for an unlabelled combinatorial class C, for which there 
are c n elements of size n, its generating function is defined as 



C(z) := 2 



c n z 



Through evaluation, all the coefficients of a generating function are smashed together, and the resulting 
probabilities take into account objects of all sizes. Thus, while you do know that the object returned will 
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be uniformly sampled among objects of the same size, the size itself is a random variable — which you have 
no direct control over. As a result, a significant aspect of Boltzmann sampling involves: rejecting objects 
which are not within the desired size interval; manipulating the generating functions so the size distribution 
is such that not too many objects need be rejected. 

1.1 Obstacles to widespread popularity 

The efficiency of this approach, combined with its mathematical appeal — in many regards Boltzmann sam- 
pling is an elegant and natural application of Analytic Combinatorics pioneered by Flajolet and Sedge- 
wick [8] — have made it a fertile topic, and many of its aspects have been developed through many papers. 
Yet somehow it seems to not have caught on outside of the immediate community of Flajolet. 

The Boltzmann model. With a Boltzmann sampler for an unlabelled combinatorial class C, for which 
there are c n elements of size n, the probability of drawing an object 7 e C is 

ItI °° _ 

P,[7] = ^ with C(z) := ^ c n z n = 2 *H 

^^ Z > n=l 7 eC 

where I7I denotes the size of object 7 and x is some control parameter to be chosen. Thus the probability of 
drawing an object of size n is 

r 7 n 1 

P Z [W=n] = ^v P s [ 7 | |7l = n] = ~ 

C{Z) C n 

while the probability of drawing an object conditioned on its size is uniform. 

The name of the method evokes the Boltzmann model of statistical physics that assigns to each possible 
state of a system probability e^ /3E /Z, where E is the energy of the state, /3 = 1/T is a constant, and 
Z is normalizing constant. It is possible that the name and analogy to physics have unduly hindered the 
popularization of the method. 

Furthermore, in truth, the distribution of the objects is a very generic distribution already known to 
probabilists as the Power Series Distribution^ and according to Nelson et al. ifTOl §2.2], the terminology is 
usually credited to Noack lfl3l in 1950. 

Evaluating generating functions near their singularity. Boltzmann samplers depend on the evaluation 
of generating function in the neighborhood of their singularity — what more in constant time arithmetic 
complexity. 

Time has shown this to be problematic for two reasons. First of all, although the generating functions are 
far from exotic owing to their combinatorial origins, it seems particularly counter-intuitive to evaluate them 
near their singularity, precisely where they are likely to behave badly. Second, because explicit generating 
functions (that is which have no closed form in terms of simple functions) are rare, it is not practical to 
evaluate them in traditional programming languages and thus up until now, Boltzmann samplers have mostly 
been confined to computer algebra systems, which have the means to manipulate complex functions. 

To a degree these issues have been addressed by Pivoteau et al. lPl4l [131 . who developed a variant of 
Newton's iteration for combinatorial systems. But this iteration, while plausibly an optimal way to evaluate 
such generating functions, is quadratic. Thus its practicality and efficiency for Boltzmann sampling has 
yet to be put to the test, as it has never been implemented in a programming language — the only known 



'The Poisson, geometric, log-series distributions are all special cases of this distributions, a fact which is put to use by Flajo- 
let et al. [7. §2] who designed the Von Neumann/Flajolet scheme to simulate power series distributions using only random bits. 



prototype to partially automatize Boltzmann sampling was introduced by Canou and Darrasse [2], but it is 
limited to sum and product (i.e., algebraic types) which only require finitely many constant be computed, 
unlike more complicated operators such as the multiset and cycle constructions [6]. 



1.2 Our contribution: an extended, practical framework 

Our idea appeals to the classical random generation concept of rejection (see for instance Devroye's chapter 
on the rejection method [3. §2]). Instead of putting so much effort into evaluating the generating functions, 
we pick some nearby point that is easier to computer. 




Figure 1. Plot associated with the combinatorial specification of binary trees where all nodes are counted, 
B = Z + Z x B 2 . The thick black curve at the bottom plots C(z) = z + zC(z) 2 , or all coordinates (z, C(z)) 
usually considered for Boltzmann sampling; the shaded area is the region verifying c > z + zc 2 , and from 
which we get coordinates (z,c) which we use in our modified model. 

This is illustrated in Figure[T] Both Boltzmann samplers and our samplers use coordinates from the sha- 
ded region. But while Boltzmann samplers limit themselves to the coordinates that belong to the thick black 
curve at the bottom of the region, we allow ourselves to pick any point within the region. Of course, this 
comes at the cost of additional rejection, but we show that this rejection is constant and that for reasonable 
choices of coordinates it is practically negligible. 

By not restricting ourselves to a fixed curve, we are given some latitude with which to pick the points. For 
instance, rational numbers yield probabilities that are much easier to simulate exactly and efficiently [IT], 
and are easier to manipulate exactly (notably when using dichotomic search). Since rational numbers are 
dense in real numbers, it will prove very useful to only use rational numbers in our samplers — and this is a 
freedom of action that is afforded us by our use of rejection. 

Finally, by not having to evaluate generating functions anymore, we avoid any problems related with 
having to compute or represent these possibly implicit functions in the absence of a computer algebra system 
or library. Thus, our idea paves the way for the development of a full standalone system providing Boltzmann 
samplers. 

The ideas presented here followed from the first author's work to extend Boltzmann samplers to infinite 
objects [ 1 ], and the second author's attempts to modify the generating functions to shape the size distribution 
of the sampled objects. 



2 Analytic Samplers 

In this section we give the main definitions for our analytic random samplers, and then show the algorithms 
associated with the basic constructions. 

2.1 Main definitions 

DEFINITION 1 . Let A be an unlabelled combinatorial class, and a n the number of objects from A that have 
size n. The ordinary generating function (OGF) associated with class A is defined equivalently by 



A{z) := J] a n z n or A{z) := J] 

n=0 aeA 



*N 



The ordinary generating function enumerates combinatorial class it is associated to. The tenet of the 
symbolic method (H is that if a combinatorial class can be symbolically specified using a set of operators 
(disjoint union, Cartesian product, sequence, multiset, etc.) from initial terminal symbols called atoms 
which have unit sizdj then this specification can be directly translated to obtain the ordinary generating 
function. 

DEFINITION 2. Let A be a symbolically defined combinatorial class which can be translated, following the 
symbolic method, to afunctional equation on A{z), the generating function associated with A, 

A = $(Z, A, X) => A{z) = (f)(z, A(z),X(z)), 

where both <£> and (/) may possibly involve other classes/generating functions which we note using vectors in 
bold (and each symbol/generating function component of the vector itself defined by their own equations). 
A pair of coordinates (z, a) is said to be analytically valid coordinates for the combinatorial class A if and 
only if they verify the inequality 

a > (f)(z, a,x). 

In general, for convenience and clarity, we will omit the vector in the notations, and any additional bound 
symbols will be implicit. 

DEFINITION 3. An analytic sampler for an unlabelled combinatorial class A is an algorithm which samples 
an object a 6 A, of size \a\, with probability 

z \<*\ A(z) 
P(2 >a ) [a] = and fails with probability P(z,o) [fl = 1 

where A(z) is the ordinary generating function associated with class A, and the analytically valid coordi- 
nates (z, a) are called the control parameter. Moreover we denote by TA(z, a) such an analytic sampler. 

Theorem 1. Let A be a combinatorial class and A(z) its generating function, and let (z, a) be analytically 
valid coordinates for A. The proportion of objects of which the generation has failed is constant, does not 
depend on the size of the finally outputted object, and is equal to A(z)/a. 



2 Think of atoms as being nodes in a tree, in a graph, elements of a permutations, balls in an urn, etc. 



Proof. This follows from the definition of the model of Analytic Samplers wherein the probability of a 
single draw failing is constant — in the sense that it does not depend on the size of the object that was being 
constructed when the sampling failed — and equal to 

tp m i A( - z) 

P(*,a)l7J = X ~—' 

and thus an object (of some random size) is drawn with the complementary probability. The number of 
failures before an actual object is drawn is then geometrically distributed with p = 1 — A(z)/a. We then 
have: 



fe=0 

A{z) 



a 



A(z) 



a ( A(z) 



A(z 

a 



1 



A(z) 

and because there is one last object generated (the one that does not fail, and after which we are done) the 
expected proportion of objects which have failed is l/(Kf za \ [#+] + 1) as stated. □ 

Note that the lower bound for a is a = A(z), and for this choice of value, the analytic sampler does not 
fair] indeed the inequality can naturally be seen as an equality involving a slack variable 5, a = <ft(z, a, x) + 
5, where 5 is the proportion of failures. In essence, if you are willing to spend the computational time needed 
to compute the generating function, then you are rewarded for your efforts by having no rejection at all. 

Example 1. At this point we will illusttate these definitions by looking at the class B of binary trees, in which all 
nodes both internal and external count towards the size of the tree. These trees can be symbolically specified as either 
a leaf (Z) or node which has two subtrees (Z x B 2 ), 

B = Z + ZxB 2 and B(z) = z + zB(z) 2 . (1) 

This functional equation can then be translated to an inequality, 

b > z + zb 2 . 

The analytically valid coordinates for B are all points that belong to the shaded region in Figure [T] 

2.2 Constructions 

In this subsection, we give the basic constructions used by our analytic samplers. We follow the notation of 
the original article [5], which we extend it to include our failure probability, 

TA : [pi] ■ Berfe) => X \ Y 

means that we first fail with probability 1—pi, then we draw a Bernoulli variable U of parameter p2, if it is 
equal U = 1 then we return X if not we return Y . We instead of a Bernoulli distribution, we have a discrete 
distribution K, we mean that we return a tuplet of K independent calls to the sampler. 

Let A, B and C be combinatorial classes. We recall, for clarity, that in this article we note ¥[a] the 
probability of drawing an object a; when we want to make explicit from which class this object is drawn, 
we note P[aei]. 



3 Because Boltzmann sampler/analytic sampler already involve some kind of rejection to target the size, we call this second type 
of rejection "failure" to avoid any confusion. 



Disjoint union. Let A = B + C, and ao ^ &o + c o- We first reject 1 — (60 + co)/ao of the objects, then 
we do a normal Boltzmann. 



TA: 



bo + c 
a 



Ber 



bn 



\bo + c 



TB I TC 



(2) 



Proof. We must show that the sampler TA returns objects a e A with the correct probability P( ZjO0 ) [a] = 
z' a l/ao (that is the probability of drawing an object from A follows the law given in Definition 2), assuming 
inductively that the generators TB and TC are correct. 
Hence: 



Wa^H^^W^B] 



CQ 

bo + c 



(2. co) 



[aeC] 



That is the probability of sampling an object from A is the probability of first not failing, (bo + cq)/oo, 
and then the probability of drawing the object using the sampler for class B or C with the correct Bernoulli 
probability. By hypothesis those two samplers return objects with correct probability, so 



\z,ao) [ a e -4] 



bo + c b z 



M 



c z 



M 



a 



fro + co b b + c cq 



M 



fl 



Note that in this proof, and the following, we do not explicitly prove the probability of failure as it is a 
straightforward consequence: sum the probability of drawing an object over all possible possible objects, 
and take the complimentary probability. □ 



Cartesian product. Let A = B x C, and ao > b ■ c . 

bo -co 



TA 



«o 



(TB ; TC) 



Proof. The proof follows the same model as the previous construction; let a = (/3, 7), 



(z,ao) 



[a e A] 



bo -cq 
a 



\ zM [(3 e B]F M [ 7 e C] 



and since the samplers for TB and TC are inductively assumed to be correct, 



\z,a )[a^A] 



a b cq 



('Hi 



('Hi 



(3) 
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Sequence. Let A = Seq (B) and a ^ 1/(1 - b ). 



TA: 



a 



Geo(b ) => (TB, 



Proof. We follow the same model as previously. Let a = (fii, ...,/3k)- 

._! k 



\z,a ) [" e A] 



1-6, 



oj 



Of) 



»[Geo(6 ) = fc]niP(,,M[A eB ] 



i=l 



(4) 



and by hypothesis 



(1 _ h ) _1 k r\fo\ 

V, ao )[« £ X] = ^ ^&o fc (l - k) n — 

fl-ftnW , |/3i|+...+|/3 fc | 

z |ft|+...+|/3 fc | 2 |a| 



a a 

a 

2.3 The pointing operator (differentiation of generating functions) 

A differentiated generating function is just a generating function like any other. It requires no special treat- 
ment, except we must consider it like a new generating function, with a new variable. We can illustrate this 
with pointed binary trees. 

Let B = Z + ZB 2 be the class of binary trees where all nodes are counted. Then, following the 
traditional definition of 9 the, pointing operator (which consists in the objects resulting from all the ways to 
distinguish one unit-size atom from objects of B), we have the @B = Z + ZB 2 + Z(@B)B2 + ZB(@B). 
The associated inequalities are 

b > z + zb 2 + 2zbb 
b^ z + zb 2 . 

We then make a sampler for QB using the previously described constructions. 

3 The substitution operator 

The substitution (or composition) operator [8, §1.6.2] can be seen as replacing every of one class with 
elements of another. As an example, the class U of unary-binary trees can either be defined traditionally as 
either a leaf (Z), a unary node (ZU) or a binary node (ZU 2 ) or can be seen as the class B of binary trees (as 
defined in Equation ([T]), where all nodes count towards the size of the tree), in which we substitute all nodes 
by a strictly positive chain of atoms Seq^ (Z). This corresponds to the specifications 

U = Z + ZU + ZU 2 or U = Bo Seq m (Z) (5) 

which respectively translate to 

U(z) = z + zU{z) + zU{z) 2 or U(z) = B ( -—] . (6) 



The substitution operator is also central in other constructors such as the cycle, the set or the multiset — 
collectively referred to as the Polya operators, because they require carefully considering symmetries and 
making sure these are not counted multiple times. 

How to substitute when there are no functions? It might seem at first glance that our model, where all 
generating functions are replaced by fixed values, by essence procludes any such substitution. Such is not 
the case. Remember that an analytic sampler for class A is parameterized by a set of coordinates (z, a); and 



indeed the first coordinate acts can be modified much like the value in which the generating function would 
be evaluated. 

Formally, let B and C be two combinatorial classes for which we are able to construct analytic samplers, 
and let A be a class defined by substitution, A = B o C. An analytic sampler for A can be built as follows: 
let (z, c) be analytically valid coordinates for C and let (c, b) (note that we reuse the same variable c as the 
first coordinate) be analytically valid coordinates for B. Then 

3.1 A simple example of unordered operator 

In this subsection, we guide the reader through the use of the substitution operator by giving two examples: 
first we show how to realize samplers for a modest operator, MSET2, for the unordered pair; then we use 
this to exhibit a sampler for non-plane binary trees, also called Otter trees. 

Unordered pair. Let B be a combinatorial class, and A = MSET2 (B) be the class containing unorde- 
red pairs of elements of B. The corresponding generating functions A{z) and B{z) verify the functional 
equation 

A{z) _ Bj^m^ (7) 

Assuming there is an analytic sampler for B, we can build an analytic sampler for A. Let (z, b) and (z 2 , b) 
both be analytically valid for B (note that the variable z must be the same in both pairs), and let (z, a) be 
analytically valid for A, that is 

b 2 + b 
a > —j-. (8) 

Using the notation we have introduced, 



TA(z,a) : 



b 2 + b 
2a 



2 



Ber( -^— = ) => {TB(z,b) ; TB(z,b)} \ TB{z 2 ,b) and duplicate. (9) 



In other terms, after making the obligatory failure test, we choose with the proper probability whether 
to create a pair of elements resulting from independent calls to TB(z, b), or whether to make one call to 
TB(z 2 , b) and duplicating the resulting object to make a pair of identical objects. 

Proof. As before, proving the validity of this algorithm involves showing that the analytic sampler TA (z, a) 
returns any object a e A with probability zH/a. We distinguish two disjoint cases. 

• Either the pair a = {Pi; P2] contains two distinct elements, /3i 4= /?2- Then this pair could only have 
been produced by two independent (and distinguished) calls to TB(z, b). Thus under this setting, 

P 2 ,a[a I 01 + H = b ~l^ ■ 7^ <Pz,b\Pl] ■ WzM + V z ,b[fo] ■ Pz, 6 [/3l]) ■ (10) 

2a b 2 + b 

By hypothesis, TB(z, b) is an analytic sampler for class B, which means it returns an object (3 e B 
with probability z\P\/b, 

h 2 + l h 2 2 *I*I + IAI Z \h\ + \M z \a\ 
W z , a [a ft + fa ] = ~^— ■ To— r r 2 = = " (11) 

2a b 2 + b a a 



• Or the pair a = {/?; /?} contains two identical objects. The pair could then have been drawn by either 
branch: from two independent calls to TB(z, b) which happen to return the same object; or from the 
call to FB(z 2 , b) which is duplicated. In this case, 

p„[« i /,- a _ *±» . (^ . r ,M + ^ ■ p«, s ra) ■ <i2> 

Again assuming the analytic sampler for B is correct, 

P 2 , a a /? = /? = — — • — -= ■ — + — -= • U=U = = . (13) 

2a \ b 2 + b\bj b 2 + b b a a 

a 

Otter trees. We've already discussed thoroughly discussed the class B of binary trees. These binary trees 
are plane, in the sense that there the children of an internal node are distinguished: there is a left node and a 
right node. We now consider the class O of Otter tree, which are binary trees that are non plane, using our 
newly defined MSET2 operator, 

O(z) 2 + Otz 2 ) 
= Z + MSet 2 {0) => 0{z) = z+ { ' { ' . (14) 

Note that we only count external nodes. This translates to the initial set of inequalities, 

O[ ] 2 + O[i] ^ 2 °[i] 2 + °[2] 

°[0] > z + 2 °[i] > z + 2 ■ 

The analytic samplers are then 

z+ (°[0] 2 + O[i])/2 



TO(z,o [0] ): 



D 



°[0] 

2 



•Ber 



z+ (o[0] 2 + O[i])/2 



Ber —^- I =» (TO(z,o [0] ) ; TO(z,o [0] )) \ TO(z 2 , 0[1] ) and duplicate 

V°[o] +°[i]/ 



Unsurprinsingly each level of recursion for the tree requires squaring the value of the current z, and also 
calculating a new ou — implemented in the existing Boltzmann model, this would have instead required 
reevaluating the generating function. As it turns out, the om must be deduced from the inequalities for 
i < 2, but can then be deduced automatically. An implementation is to follow in the full version of this 
paper. 

3.2 Multiset and other Polya construction 

The reason we have taken such care in describing a few simple applications of the substitution operator, 
is so the reader may get a feel for the way it works. But one of the advantages of this substitution ope- 
rator is in sampling from objects containing more advanced constructions such as the set, multiset, cycle, 
which involve symmetries as we've hinted at previously, and thus generating functions which require heavy 
computation. 




Indeed as an example, let A be a combinatorial class. The class M. = MSet (A) of multisets of 
elements of A has the following generating function 

M(z) -- 

Typically the infinite sum is of course truncated, to the first twenty terms for instance [6], but even then 
it requires evaluating the generating function A(z) in as many points. Instead with our method, we only 
require calculating the exponentiations z, z 2 , z 3 , etc. and finding some good enough values. 

with probability 1- (exp (sum (ai/i, i=l .. infinity) ) /tO) fail 

M <- {} 

p <- Maxlndex() 

for i = 1 to kO - 1: 
ki <- PoissonLaw (ai/i) 
for j = 1 to ki : 

<- Gamma (A) (z A i, ai) 

M <- M u { 0, 0. . . } replicated i times 
kp <- PositivePoissonLaw (ap/p) 

for j = 1 to k (kO) : 

<- Gamma (A) (z A p, ap) 

M <- M u { 0, 0. . . } replicated p times 

return M 

Proof. Let m = {ai kl , a^ 2 , . . . , a n kn } be the multiset's composition. We consider all possible ways u to 
draw m, which note ar(m). 

For instance, if m = {a, a, a, b, b} which we note m = {a 3 , 6 2 }, the ar(m) arrangements are: a|a|a|6|6, 
a|a|6|a|6, a|6|a|a|6, 6|a|o|a|6, a|a|6|6|a, a|6|a|6|a, 6|a|a|6|o, a|6|6|a|a, 6|a|6|a|a, 6|6|a|a|a; a\a 2 \b 2 ; a\b\b\a 2 , 
6|o|6|a 2 , 6|6|a|a 2 ; a|a|a|6 2 ; b 2 \a 3 ; b\b\a 3 . The arrangement a|6|fo|a 2 corresponds to having drawn p = 2 (the 
maximal index), p\ = 3, and then the objects a,b, and b, then p2 = 1 and drawn b (which is thus then 
replicated twice). The arrangement 6 2 |a 3 corresponds to having drawn p = 3, p\ = 0, p2 = and b draw 
and replicated twice, and P3 = 1 and a drawn and replicated three times. 

Thus an arrangement can be written as 

/ «(1), v(2), , v(f)\ 

{ a u(l)\ a u(2)\---\ a u(f)J 

with v(l) 5= v(2) > «(/) and for all 1 ^ i > n, T,k;u(k)=i v ^ = ki ~ 

We immediately notice that pj = \ {i, v(i) = j} \ and that {i, v(i) = j} is an interval which we will call 



';• 



We now must calculate for all arrangement u the probability with which it is obtained by the algorithm: 
• we must draw the maximal index p := P(K = p) = P(K < p) — P(K < p — 1) thus if p = 0, 
P(K = = ^ — — and p > 0, P(K = p = ?—= — . 

1 



then we must draw all the pi (following a Poisson law) for i < p which gives fT^ {ti/i) Vi exp(— ti/i) 

10 



Pi\ 



then draw p p (following a positive Poisson law) which yields {t p /p) Pp (exp(t p /p) — 1) x — - 

Pp- 



• 



then all the calls to the analytic sampler which yield 



z, 



\m\ 




rw 



z \m\ 1 

After simplification, we are left with t==-^ — Y,„, c „„,( m \ fT? -* — ; — ■ We are left with having to 

show that ZuUj —^ = 1- 

To that affect, we associate with each arrangement u, the monomial M(u) = \\ s^ 3 where 

<Hj = \{k,v(k) = j and u(k) = i}\. 

For example, for u = a\b\b\a 2 , the corresponding monomial is M(u) = si,is| \Si,2 (si,i corresponds to the 
a, «2 i to the two b and si^ to the a 2 ). 

To conclude, we must now point out that 

n -. 

n[^]exp(2K fc x fc A))= 2 M( U )[]-^ 

i=l uear(m) j - 1 

avecm = {c^ 1 , a^ 2 , ...,a^ n }. 

Then, by evaluating both sides of the equation with the Sij set to 1, we finish the proof. □ 

4 Determining analytically valid coordinates 

In this section, we finally discuss how we obtain values from the system of inequalities translated from 
the specification. The system of inequalities, because of their combinatorial origin, has a certain type of 
"convexity": this is the convexity witnessed on Figure 1, which generalizes to higher dimensions, and 
which we characterize by the following two lemmas. 

Lemma 1 (Initial segment). Let Abe a combinatorial class, and let (z, a) be analytically valid for A. Then 
for all y e [0, z], the coordinate (y, a) is also analytically valid for A. 

Sketch of proof. Let A{z) = (f>(z, A(z), X(z)) be the functional equation associated with A as stated in 
Definition 2; because <fr(z, a, x) is analytic in (0, 0), it admits a series development, 

<p(z, a,x)=Yj c i,j,k z % a?x k (15) 

i,j,k 

because of its combinatorial origin (it is built only from additions, products, exponentiations, etc. which 
preserve positivity), the coefficients Cij^ are all positive. It is now straightforward to see that 

a > (f>(z,a, x) = Yj Ci,j,k z % a j x k => Vy 6 [0, z],a > 4>{y, a, x). (16) 

i,j,k 

□ 

Lemma 2 (Vertical convexity). Let A be a combinatorial class, and let (z, a 1 ) and (z, a") both be ana- 
lytically valid for A (and assume without loss of generality that a' < a"). Then for all a e [a' , a"], the 
coordinate (y, a) is also analytically valid for A. 

11 



Sketch of proof. Consider a = a'5 + a"(l — S) with 5 6 [0, 1]. By convexity ofx^ x?, we have 

(a') j S + (a") j (l-6)> {a'5 + a" {I - 8)) j . (17) 

Thus 

a > J] c id,k z l Ua'Y 5 + (a") j (1 - <5)1 x k > ^ c m> z * ( a ' J + a "( 1 ~ S )Y x * = ^( z > a ' x )- ( 18 ) 

a 

Because of this convexity, it is now possible using a naive dichotomic search algorithm to locate the 
best value for one variable in logarithmic time. If there are k values to find, we expect that searching for 
each of them one after the other independently (as allowed by Lemma 2), we would obtain an algorithm 
in O(klogi) where t is the precision with which we want to determine the values. It remains to be seen 
whether we can improve on this naive idea. (For the moment, in our implementations, values are entered by 
hand.) 

5 Conclusion 

The recursive method was probably the first systematic method to automatically build samplers from a sym- 
bolic specification, but it required a lot of preprocessing. Boltzmann sampling was designed with the idea 
to avoid any preprocessing by using the generating function has a global source of the mixing probabilities. 
But it turns out that evaluating the generating functions could be overly precise a problem for the way it is 
used. 

In this paper, we propose to use the classical idea of rejection to relax the way that the probability 
mixing is done. In this way, we have more latitude with the values we pick for the samplers. This does 
not completely make computations disappear as we still have a system of inequalities to solve, but this is a 
relatively simpler problem to solve, and we believe it is much more flexible. One example is that we can 
limit ourselves to using rationals which is a great advantage for the simulation of the random laws. 

Furthermore, while we are only at a tentative stage, and our proposition is, for the time being a bit rough, 
many extensions seem possible, notably the developments in Boltzmann sampling which evolve changing 
the control parameter: in Boltzmann sampling this would mean reevaluating the generating function many 
times, but as we showed, in our model the substitution only involves calculating the parameter and no 
evaluation is required. 
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